%===============================================================
%=        Noise Cancellation With Adaptive Filters             =
%===============================================================
%=   The Coeficients Of The Adaptive Filter Are Calculated     =
%=                       Using The                             =
%=          RECURSIVE LEAST SQUARES (RLS) ALGORITHM            =
%=     The Algorithm Developed First By Plackett In 1950.      =
%===============================================================
%= 						        Note                              =
%=                       ------------                          =
%= In The Variable FF (Forgetting Factor) Is Stored The Value  =
%= That Specifies The Forgetting Factor For The RLS Algorithm. =
%= Typically, FF Is Between 0.98 and 1. Smaller Values Assign  =
%= Too Much Weigth To The More Recent Data, Which Leads To     =
%= Wildly Fluctuating Estimates.                               =
%===============================================================

clc
clear
%The Forgetting Factor.
FF = 0.98; 								                 
%The Filter Length.
FLength = 40; 							                 
%The Number Of Iterations.
NIterations = 1000; 					  
%The Input White Noise.
Noise = 0.7 * randn(NIterations,1);

for n=1:NIterations 
   %The Input Signal.
   InputSignal = 2*sin(pi/2*linspace(1,50,1000));
end

figure (1);
plot(InputSignal,'b');
title('The Input Signal');
ylabel('Amplitude');
xlabel('Time');
grid on

%The Signal Plus The Noise.
SignalPlusNoise = InputSignal' + Noise;

figure (2);
plot(SignalPlusNoise,'b');
title('The Signal Plus Noise');
ylabel('Amplitude');
xlabel('Time');
grid on

%========================================================
%= The Adaptive Algorithm RLS (Recursive Least Squares) =
%========================================================
%= START =
%=========

%Initialise The Filter Parameters By Setting:
INoise = 0*(1:FLength)'; 
RYY = linspace(0.01,0.01,1000); 
FCoeff = linspace(0,0,40);

%Calculation Of The First Iteration.
%Setting All The Variables To The Initial Values.
OSignal = FCoeff * INoise;
Error = Noise(1) - OSignal;
MSE(1) = Error;

Number = RYY(1);
Nominator = RYY(1) * INoise' * RYY(1) * INoise;
Denominator = FF + INoise' * RYY(1) * INoise;
RYY(2) = (Number - (Nominator / Denominator)) / FF;
FCoeff = FCoeff + (RYY(2)* INoise * Error)';

%Calculation Of The Nth Iteration.
for j = 1:NIterations
   
INoise(1) = Noise(j);                         
OSignal = FCoeff * INoise;
Error = SignalPlusNoise(j) - OSignal;
MSE(j) = Error;

Number = RYY(j);
Nominator = RYY(j) * INoise' * RYY(j) * INoise;
Denominator = FF + INoise' * RYY(j) * INoise;
RYY(j+1) = (Number - (Nominator / Denominator)) / FF;
FCoeff = FCoeff + (RYY(j+1)* INoise * Error)';

end 

%=======
%= END =
%=======================================================
%= The Adaptive Algorithm RLS (Recursive Least Squares)=
%=======================================================

figure (3);
plot(MSE,'b');
title('The Filtered Signal');
ylabel('Amplitude');
xlabel('Time');
grid on

%The Frequency Spectrums For The Unfiltered and Filtered Signals.

Time = 0:8000/(512-1):8000;
figure(4);
plot(Time,abs(fft(SignalPlusNoise,512))./512,'b');
title('Fast Fourier Transform (Frequency Domain) For The Noisy Signal');
xlabel('Frequency');
ylabel('Magnitude');
grid on

figure(5);
plot(Time,abs(fft(MSE,512))./512,'b');
title('Fast Fourier Transform (Frequency Domain) For The Filtered Signal');
xlabel('Frequency');
ylabel('Magnitude');
grid on

%The Time Domain For The Mean Square Error

FreqXAxis = 0:length(SignalPlusNoise)-1;
figure(6);
plot(FreqXAxis,(InputSignal-MSE).^2,'b');
title('The Mean Square Error');
xlabel('Time');
ylabel('Magnitude');
grid on

%The Final Filter Coeficients, The MSE And The Signal-To-Noise Ratio

disp('Estimated Filter Coefficients FCoeff =');
disp(FCoeff');

disp('Mean Squared Error (MSE) Error =');
disp(abs(Error));

SignalToNoiseRatioUnfiltered = 10*log10(sum(abs(fft(InputSignal)).^2)/sum(abs(fft(SignalPlusNoise-InputSignal')).^2));
disp('Signal To Noise Ratio For The Unfiltered Signal =');
disp(SignalToNoiseRatioUnfiltered);

SignalToNoiseRatioFiltered = 10*log10(sum(abs(fft(InputSignal)).^2)/sum(abs(fft(MSE-InputSignal)).^2));
disp('Signal To Noise Ratio For The Filtered Signal =');
disp(SignalToNoiseRatioFiltered);
